perm filename FFFT.IL[TIM,LSP] blob sn#722267 filedate 1983-07-28 generic text, type T, neo UTF8
(FILECREATED "22-FEB-83 04:12:43" {PHYLUM}<GABRIEL>FFFT.;3 3984   

      changes to:  (FNS FFFT)
		   (VARS FFFTCOMS)
		   (MACROS IEXPT)

      previous date: "22-FEB-83 03:48:03" {PHYLUM}<GABRIEL>FFFT.;1)


(* Copyright (c) 1983 by JonL)

(PRETTYCOMPRINT FFFTCOMS)

(RPAQQ FFFTCOMS ((FILES (SYSLOAD FROM LISPUSERS)
			CMLARRAY)
		 (FNS FFFT)
		 (VARS (RE (MAKEARRAY 1025 (QUOTE INITIALELEMENT)
				      0.0))
		       (IM (MAKEARRAY 1025 (QUOTE INITIALELEMENT)
				      0.0)))
		 (MACROS IEXPT)))
(FILESLOAD (SYSLOAD FROM LISPUSERS)
	   CMLARRAY)
(DEFINEQ

(FFFT
  (LAMBDA (AREAL AIMAG)                                      (* JonL "22-FEB-83 04:11")
                                                             (* Fast Fourier Transform AREAL = real part AIMAG = 
							     imaginary part)
    (PROG (AR AI PI I J K M N LE LE1 IP NV2 NM1 UR UI WR WI TR TI)
          (SETQ AR AREAL)                                    (* Initialize)
          (SETQ AI AIMAG)
          (SETQ PI 3.141593)
          (SETQ N (ARRAYDIMENSION AR 0))
          (add N -1)
          (SETQ NV2 (LRSH N 1))
          (SETQ NM1 (SUB1 N))
          (SETQ M 0)
          (SETQ I 1)
      L1  (until (NOT (ILESSP I N))
	     do                                              (* Compute M = log (N))
		(add M 1)
		(add I I)
	     finally (if (NOT (IEQP N (IEXPT 2 M)))
			 then (PRINC "Error ... array size not a power of two.")
			      (READ)
			      (RETURN (TERPRI))))
          (SETQ J 1)                                         (* ;Interchange elements)
          (SETQ I 1)                                         (* ;in bit-reversed order)
      L3  (repeatuntil (NOT (ILESSP I N))
	     do (if (ILESSP I J)
		    then (SETQ TR (\PAREF AR J))
			 (SETQ TI (\PAREF AI J))
			 (\PASET (\PAREF AR I)
				 AR J)
			 (\PASET (\PAREF AI I)
				 AI J)
			 (\PASET TR AR I)
			 (\PASET TI AI I))
		(SETQ K NV2)
		L6
		(until (NOT (ILESSP K J))
		   do (SETQ J (IDIFFERENCE J K))
		      (SETQ K (LRSH K 1)))
		(SETQ J (IPLUS J K))
		(add I 1))
          (for L to M
	     do                                              (* ;Loop thru stages)
		(SETQ LE (IEXPT 2 L))
		(SETQ LE1 (LRSH LE 1))
		(SETQ UR 1.0)
		(SETQ UI 0.0)
		(SETQ WR (COS (FQUOTIENT PI (FLOAT LE1))))
		(SETQ WI (SIN (FQUOTIENT PI (FLOAT LE1))))
		(for J to LE1
		   do                                        (* ;Loop thru butterflies)
		      (for I from J by LE until (IGREATERP I N)
			 do                                  (* ;Do a butterfly)
			    (SETQ IP (IPLUS I LE1))
			    (SETQ TR (FDIFFERENCE (FTIMES (\PAREF AR IP)
							  UR)
						  (FTIMES (\PAREF AI IP)
							  UI)))
			    (SETQ TI (FPLUS (FTIMES (\PAREF AR IP)
						    UI)
					    (FTIMES (\PAREF AI IP)
						    UR)))
			    (\PASET (FDIFFERENCE (\PAREF AR I)
						 TR)
				    AR IP)
			    (\PASET (FDIFFERENCE (\PAREF AI I)
						 TI)
				    AI IP)
			    (\PASET (FPLUS (\PAREF AR I)
					   TR)
				    AR I)
			    (\PASET (FPLUS (\PAREF AI I)
					   TI)
				    AI I))
		      (SETQ TR (FDIFFERENCE (FTIMES UR WR)
					    (FTIMES UI WI)))
		      (SETQ TI (FPLUS (FTIMES UR WI)
				      (FTIMES UI WR)))
		      (SETQ UR TR)
		      (SETQ UI TI)))
          (RETURN T))))
)

(RPAQ RE (MAKEARRAY 1025 (QUOTE INITIALELEMENT)
		    0.0))

(RPAQ IM (MAKEARRAY 1025 (QUOTE INITIALELEMENT)
		    0.0))
(DECLARE: EVAL@COMPILE 

(PUTPROPS IEXPT MACRO (X
  (PROG ((N (CAR (CONSTANTEXPRESSIONP (CAR X))))
	 (E (CADR X)))
        (RETURN (if (AND (FIXP N)
			 (POWEROFTWOP N))
		    then (if (NEQ 2 N)
			     then (SETQ E (BQUOTE (ITIMES , (SUB1 (INTEGERLENGTH N))
							  ,E))))
			 (BQUOTE (MASK.1'S , E 1))
		  else (BQUOTE (EXPT (IPLUS 0 , (CAR X))
				     (IPLUS 0 , (CADR X)))))))))
)
(PUTPROPS FFFT COPYRIGHT ("JonL" 1983))
(DECLARE: DONTCOPY
  (FILEMAP (NIL (560 3379 (FFFT 570 . 3377)))))
STOP